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Abstract.  We  introduce  an  approach  to  velocity  and  reflectivity  es¬ 
timation  based  on  optimizing  the  coherence  of  multiple  shot-gathering 
inversions  of  reflection  seismograms.  The  resulting  algorithm  appears  to 
avoid  severe  convergence  difficulties  reported  for  output  (nonlinear)  least- 
squares  inversion.  We  describe  in  detail  an  algorithm  appropriate  for 
plane-layered  acoustic  models,  using  the  convolutional  approximation  to 
the  plane-wave  (p-tau)  seismogram.  We  give  theoretical  and  numerical 
evidence  that  coherency  optimization,  as  defined  here,  yields  stable  and 
reasonably  accurate  estimates  of  both  velocity  trend  and  reflectivity,  by 
exploiting  reflection  phase  moveout  and  amplitudes  in  a  computationally 
efficient  way.  We  demonstrate  that  the  approach  may  be  applied  to  field 
data  by  extracting  velocity  and  reflectivity  estimates  from  a  Gulf  of  Mex¬ 
ico  marine  data  set.  Finally  we  explain  briefly  how  the  approach  may  be 
modified  to  determine  elastic  models  and  source  parameters  as  well  as  to 
determine  laterally  heterogeneous  models. 


1  Introduction 

Full  waveform  inversion  algorithms  for  multi-offset  seismic  reflect. on  data  have  been 
discussed  extensively  over  the  past  several  years.  For  a  small  sample  of  this  discus¬ 
sion,  see  Bamberger  et  al.  (1982),  Lines  and  Treitel  (1984),  Tarantola  (1984,  1986), 
Gauthier  et  al.  (1986),  MacAulay  (1985),  Kolb  et  al.  (1986),  Mora  (1986,  1987),  Pan 

°This  research  was  partly  supported  by  the  Office  of  Naval  Research  under  contracts  N00014- 
85-0725  and  N00014-89-J-1115,  and  by  the  National  Science  Foundation  under  Grant  DMS  8603614 
(W.W.  Symes).  Computations  were  performed  at  the  Cray  XMP  facility  of  the  Naval  Research 
Laboratory.  We  are  grateful  for  the  permission  of  Exxon  Production  Research  Company  to  publish 
the  results  based  on  the  field  data,  and  to  allow  the  participation  of  J  .J .  Carazzone. 
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et  al.  (19S8),  Pan  and  Phinney  (1989),  Cao  et  al.  (1989).  All  of  the  work  just  cited 
is  based  on  the  output  least-squares  principle  in  which  a  physical  model  of  the  sub¬ 
surface  is  adjusted  to  minimize  the  mean  square  error  between  model-predicted  and 
data  seismograms.  This  approach  does  not  require  picked  travel-times,  unlike  reflec¬ 
tion  tomography  (Bube  et  al.  (1985),  Lines  et  al.  (1987))  and  in  principle  can  extract 
seismic  velocity  estimates  as  well  as  reflection  amplitudes,  unlike  linearized  inversion 
(Cohen  and  Bleistein  (1979),  Clayton  and  Stolt  (1981),  Beylkin  (1985),  Ikelle  et  al. 
(1988),  Beylkin  and  Burridge  (1987)).  In  addition,  any  desired  level  of  physical  de¬ 
tail  may  be  built  into  the  output  least-squares  principle,  and  it  may  also  incorporate 
non-seismic  constraints  (Lines  et  al.  (1988)).  Accordingly,  some  motivations  for  the 
interest  in  model-based  full  wave-form  inversion  are  the  possibilities  of  more-or-less 
automatic  determination  of  complicated  seismic  velocity  structures,  and  of  rational 
extraction  of  reflection  characteristics  directly  indicating  the  presence  of  hydrocarbon 
deposits. 

In  practice  it  has  proven  difficult  to  realize  the  promise  of  output  least-squares 
inversion,  even  for  synthetic  data  sets.  The  hypersensitivity  of  the  reflection  seis¬ 
mogram  to  slowly  varying  velocity  components  results  in  severe  non-convexity  of  the 
mean-square-error,  and  consequent  slow,  or  no,  convergence  of  the  (necessarily  it¬ 
erative)  solution  algorithms.  Despite  some  progress  for  layered  models  (Kolb  et  al. 
(19S6),  Canadas  and  Kolb  (1986),  MacAulay  (1986),  Cao  et  al.  (19S9)),  in  general 
quite  accurate  initial  estimates  of  velocity  are  required,  else  the  amplitude  preserving 
depth-migration  function  is  disabled  as  well  (Gauthier  et  al.  (1986),  Mora  (1987), 
Pan  and  Phinney  (1989),  Ikelle  et  al.  (1988);  Tarantola  (19S6),  Santosa  and  Symes 
(1989),  Spratt  (1987)).  That  is,  for  essentially  mathematical  reasons,  least-squares 
inversion  generally  fails  to  give  reliable  velocity  and  reflectivity  estimates. 

The  purpose  of  this  paper  is  to  introduce  a  variant  of  output-least-squares  inver¬ 
sion  which  appears  to  have  better  mathematical  properties  than  does  the  straightfor¬ 
ward  version.  This  variant  is  based  on  optimizing  the  coherence  of  multiple  inversions, 
so  we  call  it  the  coherency  method.  The  essential  idea  is  very  old,  and  is  operative 
in  ordinary  (NMO-based)  velocity  analysis  and  in  more  speculative  iterative  before¬ 
stack  migration  schemes.  Our  specific  quantification  of  coherence  seems  to  be  novel, 
and  leads  to  a  large  class  of  algorithms,  which  we  begin  to  explore  in  this  and  related 
papers. 

The  simple  version  of  the  coherency  method  used  in  this  paper  is  based  on  a 
model  of  the  plane-wave  (p-tau)  seismogram.  Two  hypotheses  explain  the  essence  of 
the  model,  and  also  hint  at  its  limitations. 

i)  Each  (plane-wave)  trace  is  the  convolution  of  a  primary  reflectivity  with  a  source 
wavelet  (plus  noise). 

ii)  The  reflectivities  for  various  slownesses  (or  incidence  angles)  are  related  in  the 
manner  appropriate  for  a  layered  acoustic  model. 

The  first  assumption  is  that  each  trace  can  be  explained  by  the  venerable  convo¬ 
lutional  model.  Thus  multiple  reflections  are  neglected.  The  second  hypothesis  makes 
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explicit  the  layered  medium  assumption,  already  implicit  in  (i),  and  also  restricts  the 
change  in  reflectivity  amplitude  with  angle  to  that  appropriate  to  a  layered  fluid. 

We  shall  establish  three  major  points  in  this  paper:  First,  a  variational  principle 
extracted  from  the  above  hypotheses  (Section  3)  is  convex  over  a  large  region  in  model 
space.  In  particular,  when  coupled  with  suitable  numerical  optimization  software, 
it  yields  the  only  waveform  inversion  algorithm  known  to  the  authors  capable  of 
reliable  estimation  of  velocity  (trends,  i.e.  long  wavelengths)  starting  from  grossly 
incorrect  initial  estimates.  We  illustrate  this  feature  via  a  synthetic  example  (Section 
4).  Mathematically  rigorous  arguments  leading  to  similar  conclusions  are  presented 
in  companion  papers,  referenced  below. 

Second,  the  coherency  method  is  practical ,  in  the  sense  that  it  can  be  applied  to 
field  data  sets  and  —  at  least  sometimes  —  produces  reasonable  results.  We  process 
a  plane-wave  gather  derived  from  a  marine  survey  in  the  Gulf  of  Mexico,  using  the 
implementation  of  the  coherency  method  mentioned  above  (Section  5).  This  data 
set  was  selected  and  processed  to  conform  as  closely  as  possible  to  the  hypotheses 
underlying  the  current  implementation  (i.e.  (i)  and  (ii)  above).  The  velocities  and 
reflectivities  obtained  compare  favorably  to  those  evident  in  a  sonic  log  from  a  well 
near  the  line. 

Third,  the  coherency  method  may  be  formulated  for  laterally  heterogeneous  mod¬ 
els,  for  variable-density  or  elastic  (rather  than  acoustic)  models,  and  using  the  fully 
nonlinear  model-data  relationship  rather  than  the  convolutional  model  or  analogous 
higher-dimensional  linearized  approximations.  In  particular,  the  coherency  method 
is  no  more  restricted  to  analysis  of  layered  systems  than  is  output-least-squares  in¬ 
version.  This  point  is  of  critical  importance,  as  the  central  issue  in  velocity  analysis 
is  reliable  and  accurate  estimation  of  laterally  heterogeneous  velocities.  We  describe 
explicitly  a  formulation  based  on  shot  gathers,  designed  to  extract  laterally  hetero¬ 
geneous  velocities  and  reflectivities  in  the  acoustic  model  (Section  6). 


2  Output  Least-squares  Inversion 

The  discussion  in  this  and  the  following  three  sections  is  based  on  the  convolutional 
approximation  to  the  plane-wave  seismogram  over  a  plane-layered  acoustic  earth 
model  with  constant  density.  Thus  the  only  mechanical  variation  allowed  in  this 
model  is  the  compressional  velocity  profile  c(z).  While  this  model  underlies  a  great 
deal  of  seismic  data  processing,  it  is  clearly  a  gross  simplification  of  seismic  wave  me¬ 
chanics.  Some  of  the  concepts  discussed  below  are  extended  to  more  realistic  models 
in  Section  6.  The  mathematical  essence  of  velocity  inversion  is  already  present  in  the 
very  simple  plane-layer  context,  however,  so  we  first  develop  the  solution  of  the  prob¬ 
lem  in  that  context.  For  discussion  of  plane-wave  seismograms,  and  their  production 
from  synthetic  point-source  data  and  field  reflection  data,  see  Gutowski  et  al.  (1982) 
and  Brysk  (1986). 

Denote  the  plane- wave  (“p-T”)  seismogram  corresponding  to  a  velocity  profile 
c(z )  by  S[c].  The  arguments  in  this  paper  are  based  on  the  well-known  convolutional 
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approximation,  which  is  reasonably  accurate  when  c  may  be  re-written  as  a  sum 

C&Cs  +  Cr 


where 


ca(z )  is  a  slowly  varying  background  velocity  model 

cr(z)  is  a  rapidly  varying  “ reflector  sequence having  locally  zero  mean 
on  the  length  scale  of  significant  change  in  the  background  velocity. 


Thus  the  (two-way)  travel-time  to  depth  z  of  a  precritical  plane-wave  at  (horizontal) 
slowness  p  is  determined  with  a  small  error  by  the  background  velocity  c4(z)  according 

r(z,p)  =  2 =2 /o  i  . 
where  vs  is  the  vertical  (plane-wave)  velocity  at  slowness  p: 


v,(z,p)  = 


^-pI) 


-1/2 


=  c4(z)A(z,p) 


A  (z,p)  =  (1  -c2a(z)p2)  1/2 . 

The  convolutional  approximation  to  the  plane- wave  seismogram  with  (possibly  directionally- 
dependent,  known)  source  wavelet  f{t,p)  is  then 


S[c„CT](t,p)  =  / 


Jq  dt'f(t-t',p)^(t',p) 


(1) 


where  the  “reflectivity”  r(t,p)  is  given  by 


r(t,p)  = 


Vr  {({t,p),p) 


by  means  of  the  inverse  two-way  travel-time  function  defined  implicitly  by 

Kb.p)  1 

t  =  2  / 

Jo  vs 

and  the  vertical  velocity  perturbation 


(2) 


Vr  =  Cr  ■  A3. 


Note  that  reflectivity  conventionally  means  dr/dt\  we  shall  confuse  r  and  its  t- 
derivative,  calling  both  “reflectivity”  as  convenient.  Note  also  that  all  surface-related 
phenomena,  as  well  as  source  and  receiver  characteristics,  are  subsumed  in  the  source 
wavelet  f(t,p).  A  detailed  derivation  of  this  model,  appropriate  to  linear  acoustics, 
may  be  found  in  Santosa  and  Symes  (1988). 
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Now  S  is  clearly  linear  in  cr,  but  quite  nonlinear  in  c3.  In  fact,  a  change  in  c3  will 
typically  result  in  a  change  in  the  “phase”  £,  and  thus  in  a  shift  in  the  high-frequency 
components  of  S,  which  in  turn  derive  from  the  high-frequency  components  in 
Since  such  a  phase-shift  may  have  a  drastic  effect  on  components  of  the  appropriate 
(high)  frequency.  Since  must  have  a  great  deal  of  high-frequency  content  in  order 
to  model  the  dense  distribution  of  reflectors  in  the  typical  sedimentary  column,  one 
expects  S  to  be  extremely  sensitive  to  changes  in  the  background  velocity  c3. 

This  oversensitivity  to  background  errors  shows  quite  clearly  in  the  expression  for 
the  perturbation  8S  in  the  seismogram,  due  to  a  perturbation  8c„  in  the  background 
velocity  (holding  cr  fixed,  and  assuming  6c,(0)  =  0): 

=  /  *  |  m,p).p) 

where 

*C(*iP)  :=  ^C(t(2,p)>p)  = 

v»(z,p)  dz'v3(z'  ,p)cf3(z')8cs(z') 

is  the  phase  perturbation  corresponding  to  8c3 ,  referred  to  depth/slowness  coordi¬ 
nates. 

Thus  the  second  derivative  of  Cr  (or,  of  r)  appears  in  the  perturbations  6S  associ¬ 
ated  with  a  background  velocity  perturbation  8c3.  On  the  other  hand,  a  perturbation 
8ct  in  cr  simply  results  in 

8S  ~  S[c5,£cv] 

as  S  is  linear  in  the  derivative  of  <v.  Thus  8S  is  linear  in  the  derivative  of  8cT. 

Again,  rv  must  be  highly  oscillatory  to  model  the  typical  reflector  distribution, 
so  the  second  derivative  of  Cj  is  typically  much  bigger,  in  any  reasonable  sense,  than 
is  the  first  derivative  of  cv  itself.  Thus  perturbation  of  the  background  velocity  c3 
has  a  much  larger  effect  on  5  than  does  perturbation  of  <v.  in  the  language  of  linear 
algebra,  the  linear  perturbation  map  (“Jacobian”) 

8c3 ,  8cr  i— >  8S 


is  ill-conditioned. 

Even  worse,  a  straighforward  extension  of  this  reasoning  shows  that  the  difference 
between  the  perturbed  seismogram 

S[c3  +  <5cs,Cr  +  8ct) 

and  its  linear  prediction 

5[cs,Cr]  +  8S 

can  easily  be  on  the  order  of  S  itself,  even  for  quite  small  8cs. 

To  summarize:  under  realistic  assumptions  on  c3  and  cr, 

(i)  the  Jacobian  8S  is  ill-conditioned; 
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Figure  1:  Background  velocity  and  reflector  series  for  synthetic  examples. 

(ii)  S  is  poorly  approximated  by  its  linearization. 

Consequently,  the  output  least-squares  objective  function 

Jls[cs,  ct]  Sdata\  :=  J  J  dt\S[cs,  cr]  —  S  fatal 

is  highly  non-convex,  with  rapidly  changing  gradient. 

To  illustrate  these  features  of  the  mean-square  error  function,  we  constructed 
the  convolutional  plane-wave  (reference)  seismogram  (Figure  2)  corresponding  to  a 
(reference)  model  (velocity  profile)  of  the  form:  smooth  background  (cs)  plus  rough 
reflector  series  (cv)  (Figure  1). 

Remark.  We  have  elected  to  calibrate  all  velocities  in  units  of  surface  velocity, 
so  that  cs(0)  =  1  always,  and  correspondingly  to  measure  depth  ( z )  in  seconds  of 
one-way  time  at  surface  velocity.  This  normalization  has  a  salubrious  effect  on  the 
scaling  of  various  numerical  computations;  physical  units  can  easily  be  recovered  when 
desired. 

We  added  multiples  of  a  trend  perturbation  Scs(z)  (Figure  3)  to  c,(z).  We  pro¬ 
duced  for  each  value  h  in  the  range  0  <  h  <  1.2  the  seismogram  corresponding  to  the 
perturbed  background  velocity 

cs(z)  +  (h-  l)6c,(z) 

with  in  all  cases  the  same  reflector  series  cr(z).  Note  that  these  are  finite  perturba¬ 
tions,  and  that  for  h  =  0  the  perturbed  velocity  profile  is  constant  (independent  of 
z).  Next  we  subtract  the  seismogram  for  h  =  1  from  the  seismogram  computed  for 
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Figure  2:  Convolutional  plane-wave  seismogram  corresponding  to  the  model  of  Fig. 
1. 

each  value  of  h ,  and  form  the  mean-square  of  the  resulting  difference  (i.e.  a  suitably 
discretized  version  of): 

[  f  dt  dp\S[cs  +  (h  -  l)6ca,cr](t,p)  -  S[ca,Cr](t,p)\2 

Jo  Jo 

This  quantity  is  plotted  versus  the  parameter  h  in  Figure  4.  This  figure  represents  a 
sampling  of  the  mean-square  error  —  the  objective  function  of  output  least  squares 
along  a  line  segment  in  model  space  connecting  the  homogeneous  background  (h  =  0) 
with  the  reference  background  ( h  =  1),  with  the  reference  seismogram  playing  the 
role  of  the  data  gather. 

The  features  mentioned  above  are  evident  in  Figure  4:  the  mean  square  error 
changes  rapidly  near  the  “solution”  ( h  =  1),  and  possesses  many  local  minima  and 
maxima.  Near  the  homogeneous  background  velocity  ( h  =  0),  which  might  be  se¬ 
lected  as  a  natural  “initial  guess,”  the  local  trend  of  the  mean  square  error  bears  no 
particular  relation  to  the  global  trend.  Consequently,  relatives  of  Newton  s  method 
have  little  prospect  of  success  in  finding  the  global  minimum.  Examples  of  the  failure 
of  Newton-type  methods  started  at  a  poor  (but  reasonable)  initial  guess  are  presented 
in  Kolb  et  al.  (1986),  Gauthier  et  al.  (1986),  Santosa  and  Symes  (1989).  Remedies 
for  this  delinquent  behavior  have  been  suggested  by  Kolb  et  al.  (1986),  Cao  et  al. 
(1989),  MacAulay  (1986).  None  of  these  remedies  change  the  basic  character  of  the 
objective  function:  highly  nonconvex,  with  rapidly  varying  gradient.  Since  methods 
based  on  local  behaviour  of  the  objective,  exemplified  by  Newton’s  method,  are  vir¬ 
tually  mandated  by  the  size  of  the  seismic  inverse  problem,  especially  for  laterally 
heterogeneous  models,  the  outlook  for  a  reliable  implementation  of  the  output  least 
squares  method  would  appear  dim. 
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Figure  3:  Velocity  trend  perturbation. 


MEAN-SQUARE  ERROR:  CONST  -  TEST.  MOD 


Figure  4:  Scan  of  mean-square  error  over  line  segment  connecting  constant  back¬ 
ground  velocity  with  that  displayed  in  Fig.  1. 
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3  The  Coherency  Method 


The  hysterical  behaviour  of  the  mean-square  error  function,  as  described  in  the  last 
section,  is  caused  by  the  interaction  between  the  background  velocity  cs  and  the 
reflectivity  r  through  the  definition  (2).  A  possible  way  out  of  this  dilemma  is  to 
decouple  c,  and  r  by  treating  r,  rather  than  tv,  as  the  “other”  component  of  the 
model:  thus  define 

5[cs,r](t,p)  =  /*  ^(<,p). 

Certainly  if  r  and  cr  are  related  by  (2)  then 

5[ca,cr]  =  5[c„r]. 

In  fact,  apart  from  the  surface  normalization,  the  background  velocity  c„  enters  the 
definition  of  S  only  implicitly,  through  the  condition  (2).  If  we  are  to  use  r  as  one  of 
the  independent  variables,  instead  of  cr,  we  must  develop  a  condition,  phrased  only 
in  terms  of  cs  and  r,  which  guarantees  that  (2)  holds.  Fortunately,  this  is  rather  easy: 
from  the  useful  perturbational  identity 

tV  _  C, 

“  <?a 

it  follows  that,  if  (2)  is  satisfied,  then 

c~z{z)cr(z)  =  u4_2(2,p)r(r(z,p),p).  (4) 

The  notable  quality  of  (4)  is  that  the  left-hand  side  is  independent  of  p.  Thus  differ¬ 
entiation  with  respect  to  p  eliminates  cr  altogether: 

0  =  ^K2(^pMt(z,p),p)].  (5) 

It  is  easy  to  reverse  this  reasoning.  Thus: 


(3) 


A  section  r(t,p)  is  the  reflectivity  corresponding  to  a  reflection  series 
Cr(z)  if  and  only  if  (5)  is  satisfied,  in  which  case  cr  is  given  in  terms  of 
r(t,p)  and  c,(z)  by  (4). 


In  formulating  the  constraint  given  by  (5),  it  is  advantageous  to  return  to  (t,p)  coor¬ 
dinates  (the  reason  will  become  apparent  below).  Thus  define  the  quantity  W[c3,  r] 
(the  incoherency)  by 


W/[ci,r](f,p) 


(z,p)r(r(z,p),p) 


z=C  (t,p) 


Then 


if  and  only  if 


f^[ca5cr]  —  S data 

■^[c*,r]  =  $data 
and  W[cs,r]  =  0 
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with  cT  and  r  related  by  (2)  and  (4). 

Of  course,  the  inverse  problem  is  overdetermined,  so  we  replace  the  exact  fit  of 
predicted  to  measured  seismograms  by  minimization  of  the  mean-square  error: 

min  J  J  |S[c„r]  -  Sdata\2dtdp 

subject  to  W[c,,  r]  =  0  (6) 

The  constrained  least-squares  problem  (6)  has  exactly  the  same  solutions  as  the 
output-least-squares  problem  discussed  in  the  preceding  section.  Since  you  can  t 
get  something  for  nothing,  something  must  be  wrong  with  (6).  The  nasty  features  of 
(6)  are  discussed  in  the  companion  papers  Symes  (1988  a,b,c).  Essentially  (6)  does 
not  have  the  properties  which  ensure  that  constrained  optimization  problems  have 
efficiently  computable  solutions,  depending  stably  on  their  data. 

To  obtain  a  problem  more  amenable  to  numerical  solution,  we  “relax”  (6)  by 
making  the  constraint  W  =  0  soft: 

min  ja[cs,  r,  Sdata]  (7) 

where 

dc[ct,r,Sdato\  =  /  /  dpdt{\S[cs,r]  -  Sdata\2  +  v2\W[c„  r]|2} 

We  call  (7)  the  coherency  optimization  problem. 

From  the  equivalence  above,  we  see  that  if  Sda^a  is  consistent ,  i.e.  Sda^a  — 
5[c,,cr],  then  the  problem  (7)  has  [c,,r]  amongst  its  solutions  for  which  r  and  cr  are 
related  by  (4).  That  is,  for  consistent  data,  (7)  has  “the  same  solutions”  as  the  output 
least-squares  problem  (6).  We  have  shown,  however,  that  (7)  is  far  better  suited  to 
numerical  computation  by  Newton’s  method  and  its  relatives,  and  that  the  solution 
of  (7)  is  stable ,  i.e.  “degrades  gracefully”  in  the  presence  of  data  error,  for  suitable 
choices  of  the  penalty  parameter  a. 

These  conclusions  hold  provided  that  Sda^a  is  near  (in  the  mean-square  sense) 
some  “exact”  or  consistent  data  S[c,,cr],  and  provided  that  c,,cT  satisfy  certain  con¬ 
ditions.  “Physical,”  or  poetic,  statements  of  these  conditions  are: 

(i)  cT  should  be  “ rough ”,  i.e.  contain  significant  variation  (reflectors),  on  a  length 
scale  dictated  both  by  the  wavelet  (/)  passband  and  by  the  smoothness  (char¬ 
acteristic  length)  of  cs; 

(ii)  The  range  of  slownesses  p  available  in  the  data  (“aperture”)  should  be  suffi¬ 
ciently  large  relative  to  both  the  degree  of  roughness  mentioned  in  (i)  and  the 
amount  of  data  error,  so  that  the  moveout  of  reflections  clearly  discriminates 
the  velocities. 

It  is  easy  to  see  why  these  conclusions  might  hold,  and  at  the  same  time  to 
identify  the  close  link  between  coherency  optimization  and  velocity  analysis ,  as  widely 
practiced  in  the  exploration  industry.  In  fact,  to  make  the  first  term  in  the  definition 
of  Ja  small,  i.e. 

J  J  dp  dt  |5[c„  r]  —  Sdata\ 
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requires  simply  deconvolving  the  data,  as  is  evident  from  the  definition  of  5,  given  near 
the  beginning  of  the  section.  Having  removed  the  influence  of  the  source  signature  /, 
the  second  term 

J  J  dp  dt\W[ca,r]\2 

gauges  the  extent  to  which  cs  reproduces  the  moveout  inherent  in  the  data.  In  fact, 
the  substitution  of  r(z,p)  for  t  in  (5)  is  a  precise  (“ray  trace”)  normal  moveout  correc¬ 
tion;  vanishing  of  the  derivative  with  respect  to  p  detects  whether  the  section  is  flat 
after  the  normal  moveout  correction.  Thus  minimizing  the  second  term  flattens  the 
tau-p  reflectivity  section  in  just  the  way  the  usual  constant-velocity  normal-moveout 
correction  flattens  events  on  a  CMP  section. 

The  coherency  method  differs  from  ordinary  velocity  analysis  in  three  crucial 
ways,  which  are  responsible  for  its  feasibility  as  an  automatic  velocity  estimator. 
First,  since  (ray-theoretic)  travel  time  is  used  to  construct  the  moveout  correction, 
rather  than  an  RMS  velocity,  the  entire  tau-p  section  is  flattened,  rather  than  small 
windows  about  certain  events.  Second,  the  reflectivity  section  (i.e.  part  of  the  model) 
is  flattened,  rather  than  data  itself.  The  derivative  figuring  in  the  definition  of  W 
guarantees  that  incoherent  error  in  the  data  will  not  be  fit  by  r,  since  such  errors 
would  greatly  increase  the  size  of  W.  For  the  same  reason,  it  is  necessary  to  apply 
W  to  the  reflectivity  model  r  rather  than  directly  to  the  data  S^a^a.  Finally,  an 
obvious  alternative  measure  of  flatness  after  NMO  correction  is  the  semblance ,  (a 
version  of)  which  would  be  generated  by  integrating  in  p  in  equation  (5),  instead  of 
differentiating,  followed  by  normalization  by  the  RMS  power.  Such  velocity  power 
spectra  are  commonly  used  detectors  of  event-flattening.  The  most  cited  reference  is 
Taner  and  Koehler  (1969).  These  integrated  measures  of  coherence  are  very  robust,  so 
can  be  applied  directly  to  the  data.  On  the  other  hand  the  effectiveness  of  semblance 
analysis  depends  on  the  non-convexity  of  the  semblance,  i.e.  on  the  sharpness  of  the 
peaks,  and  in  fact  semblance  analysis  is  essentially  equivalent  to  output-least-squares 
estimation  for  a  simplified  model  (Santosa  and  Symes,  (1989),  App.  A).  Therefore,  an 
integrated  version  of  (5)  is  unsuitable  for  automatic  velocity  estimation  by  gradient 
methods  in  the  same  way  as  is  output-least-squares. 

The  role  of  conditions  (i)  and  (ii)  above  is  clear  from  the  heuristics  of  velocity 
analysis  as  well.  The  reflector  series  Cr,  hence  each  trace  in  the  reflectivity  r,  must  be 
rough  or  oscillatory  in  the  passband  of  the  source  signature  /  in  order  that  the  tau-p 
section  contain  events  of  significant  power.  These  events  must  be  dense  enough  and 
extend  over  a  large  enough  range  of  slowness  that  the  moveout  constrains  the  velocity 
effectively  on  its  characteristic  length  scale.  In  effect,  the  typical  event  density,  seismic 
velocity  range,  and  cable  length  determine  the  resolution  limit  of  velocities.  In  the 
next  section  this  resolution  limit  will  be  illustrated  graphically. 

The  final  change  back  to  travel  time  in  the  definition  of  W  is  necessary  in  order 
that  W  define  a  smooth  function  on  appropriate  function  spaces.  Other  technical 
aspects  of  the  coherency  method,  as  well  as  precise  mathematical  statements  and 
proofs  of  our  results,  may  be  found  in  the  companion  papers  Symes  [1988a,  b,  c]. 
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H  (H  =  1:  TARGET) 

Figure  5:  Scan  of  incoherency  W  over  the  same  segment  of  velocity  models  as  used 
to  produce  Fig.  4. 

4  Numerical  Experiment:  Synthetic  Data 

Implementation  of  the  coherency  method  in  a  FORTRAN  code  (of  approximately 
3,500  executable  lines)  involved  a  large  number  of  choices  concerning  discretization, 
approximation,  and  optimization.  A  detailed  description  of  the  code  appears  in  the 
technical  report  (Symes  1988a).  For  reasons  of  space,  we  will  describe  here  only  those 
features  which  seem  to  us  most  important. 

First,  we  repeat  the  “scan”  experiment  reported  in  Section  2,  this  time  displaying 
the  coherency  optimization  functional 

Mcs,r,Sdata]  =  J  J  dp  dt  {\S[c.,r]  -  S fata?  +  <r2\W[ct,r]\2}  . 

The  integration  is  over  the  (tau,  p)  rectangle  displayed  in  Figure  2,  and  the  integral 
is  discretized  using  the  trapezoidal  rule.  The  various  derivatives  appearing  in  the 
definitions  of  S  and  W  are  approximated  by  centered  differences.  The  parameter  a 
is  set  =  1. 

The  sampling  of  over  the  line  in  model  space  (defined  in  Figures  1,  3)  is 
displayed  in  Figure  5.  The  contrast  with  Figure  4  is  dramatic,  and  provides  some 
evidence  that  the  theoretical  convexity  results  explained  in  Section  3  are  reflected  in 
the  actual  behaviour  of  J„. 

Next  we  attempt  to  recover  the  model  of  Figure  1  from  the  data  of  Figure  2,  by 
numerical  minimization  iof  Ja.  This  process  requires  the  definition  of  trial  spaces 
for  c,(z)  and  r(t,p).  For  cs(z)  we  used  a  space  of  exponentiated  integrated  cubic 


12 


6-splines: 


c,(z)  =  exp 


where 


i>j{z)  =  [N 


(*- 


zi  = 


JV 


j  =  0,...iV 


and  ^  is  a  model  cubic  6-spline  (see  DeBoor  (1978)).  This  choice  of  representation, 
parameterized  by  the  coefficients  xj ,  j  =  2, . . .  ,N—2,  ensures  positivity  of  the  velocity 
and  gives  direct  control  over  the  amount  of  smoothness  in  the  model.  See  Figure  6 
for  illustration.  (Recall  that  velocities  and  depths  are  normalized  against  the  surface 
velocity.) 


Figure  6a:  Cubic  b-splines. 
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Figure  6b:  Integrated  cubic  b-splines. 

For  r  we  used  an  array  of  samples: 

rik  «  r(iAt,pk),  t  =  0, fc  =  l,...,np 

where  At  =  .004  sec.  In  principle,  the  slowness  samples  pk  are  allowed  arbitrary  spac¬ 
ing.  In  the  synthetics  reported  in  this  section,  we  used  100  equally-spaced  slownesses 
ranging  from  0.15  to  0.35.  (Recall  that  slownesses  are  also  normalized  against  the 
surface  velocity.) 

It  is  inherent  in  our  approach  that  differentiation  of  r  with  respect  to  t  and  p 
must  be  bounded  operators,  i.e.  r  derivatives  must  be  measured  in  such  a  way  as  to 
have  essentially  the  same  size  as  r.  This  equilibration  may  be  accomplished  by  using 
discrete  analogues  of  the  Sobolev  norms  to  measure  the  size  of  r,  etc.  For  details  see 
Symes  (1988a). 

Correct  relative  weights  must  also  be  established  for  c,  and  r,  as  these  have  entirely 
different  units.  We  have  determined  the  weights  up  to  now  by  trial  and  error.  An 
approach  like  that  described  by  Kennett  et  al.  (1988)  may  circumvent  this  difficulty, 
and  is  under  study. 

We  coupled  subroutines  for  J0  and  its  derivatives  with  respect  to  the  velocity  pa¬ 
rameters  Xj  and  the  reflectivity  parameters  to  an  optimization  routine  of  truncated 
Gauss-Newton  type,  described  in  Santosa  and  Symes  (1989,  Ch.  10).  This  algorithm 
interpolates  between  steepest  descent  iteration  and  the  Gauss-Newton  iteration  in 
a  controllable  way.  In  general,  steepest  descent  iteration  is  more  reliable  far  from  a 
minimum,  whereas  Gauss-Newton  iteration  is  more  rapidly  convergent  near  a  solution 
(see  e.g.  Dennis  and  Schnabel  (1983)).  Our  algorithm  in  effect  switches  continuously 
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between  the  two,  and  achieves  assured  convergence  to  a  local  minimizer  with  some 
efficiency,  at  least  in  principle. 

Other  reasonable  choices  of  optimizer  are  certainly  possible.  We  regard  the  op¬ 
timal  choice  of  strategy  as  open,  at  present,  partly  for  the  reasons  to  be  discussed 
below. 

The  velocity  function  in  Figure  1  belongs  to  the  model  space  described  above 
with  N  =  16,  z„iaX  =  1.25  sec.  We  applied  the  algorithm  just  described  to  the  data 
of  Figure  2,  and  after  10  Gauss-Newton  steps  achieved  the  estimate  of  cs  depicted  in 
Figure  7.  The  initial  estimate  of  ca  was  c,(z)  =  1.  Rather  than  display  the  resulting 
reflectivity,  we  give  the  stack 


<%*) 


^ -  /Pm1'  dp  vJ2(z,p)r(T(z,p),p) 

iPmax  Pmin  ■'Pmin 


which  should  reproduce  Cr  if  r  is  exactly  correct.  As  Figure  8a  shows,  the  major 
features  of  cT  are  recovered,  but  with  some  loss  of  energy.  This  inaccuracy  occurs  at 
least  in  part  because  of  the  relative-scale  problem  mentioned  above.  The  estimate  for 
r  can  be  “cleaned  up”  somewhat  by  holding  c,  fixed  at  Figure  7  and  taking  a  number 
of  additional  steps  in  r  alone.  The  stack  resulting  from  this  process  is  displayed  in 
Figure  8b. 


DEPTH  (ONE-WAY  SEC.) 

Figure  7a:  Estimated  velocity.  10  Gauss/Newton/trust  region  steps;  constant  initial 
velocity,  zero  initial  reflectivity.  Shown  are  the  target  (solid  line)  and  iterates  1,  2,  3, 
and  10.  Note  that  the  general  trend  is  achieved  in  three  steps. 
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Figure  7b:  Velocity  components  of  gradient  at  iterations  1,  2,  3,  and  10. 


Figure  8a:  The  black  line  is  the  target  reflection  series  (Figure  1).  The  red  line  is 
stack  from  reflectivity  estimate  after  10  Gauss/Newton  steps. 
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Figure  8b:  The  black  line  is  the  target  reflection  series  (Figure  1).  The  red  line  is 
stack  from  reflectivity  estimate  after  30  additional  conjugate  residual  steps,  holding 
the  velocity  estimate  fixed  at  the  result  displayed  in  Fig.  7. 
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Figure  9:  NMO-corrected  final  reflectivity  estimate. 
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Another  gauge  of  the  extent  to  which  the  inversion  is  successful  is  eyeball  assess¬ 
ment  of  the  extent  to  which  the  NMO  correction  for  the  estimated  velocity  function 
flattens  the  estimated  reflectivity  section.  In  Figure  9  we  display  the  NMO-corrected 
before-stack  reflectivity 

c33{z)vJ2{z,p)r(T(z,p),p)  . 
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Figure  10:  Data  with  50%  RMS  pseudorandom  passband  noise. 


A  crude  measure  of  the  roboustness  of  this  process  is  its  sensitivity  to  additive 
noise  in  the  data  section  We  a<^ed  50%)  RMS  relative  pseudo  random  noise, 

filtered  by  the  source  wavelet,  to  the  data  of  Figure  2;  the  result  is  given  in  Figure 
10.  The  estimated  velocity  and  reflectivity  appear  in  Figures  11  and  12.  Note  the 
role  of  the  coherency  constraint  in  keeping  the  reflectivity  “clean”:  any  incoherent 
noise  drastically  increases  W,  and  is  therefore  not  modeled.  In  effect,  coherency 
optimization  implements  a  moveout-adapted  dip  filter. 


5  Numerical  Experiments:  Field  Data 

To  make  a  preliminary  assessment  of  the  practicality  of  the  method  described  in  the 
preceding  paragraphs,  we  applied  it  to  a  tau-p  data  set  produced  from  a  marine  survey 
conducted  in  the  Gulf  of  Mexico. 

The  survey  consisted  of  511  shot  profiles  (shot  interval  22.5  meters)  collected  with 
a  streamer  containing  301  hydrophone  groups.  The  group  interval  was  15  meters  with 
a  minimum  separation  of  148  meters.  Each  group  contained  17  equally  spaced  and 
equally  weighted  phones.  The  data  were  recorded  without  a  low-cut  filter;  a  110  Hz. 
high-cut  filter  was  applied;  the  sampling  rate  was  .002  sec.;  total  record  length  was 
5  sec.  A  portion  of  a  50-fold  near-trace  common  midpoint  stack  is  shown  in  Fig.  13. 
This  area  of  the  Gulf  contains  a  strong  gas-sand  related  direct  hydrocarbon  indicator 
(DHI)  located  near  2.3  sec.;  it  is  readily  apparent  in  the  stack. 

The  data  were  prepared  by  carrying  out  shot-gather  based  wavelet  estimation  and 
deconvolution  designed  to  remove  a  mild  low-  frequency  linear  noise  observed  in  the 
data.  An  isotropic  two-lobed  minimum  phase  gaussian  wavelet  was  employed.  The 
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TIME  (SEC.) 


TARGET,  0%  VS.  50%  RMS  NOISE 


Figure  11:  a.  Black  line:  target  velocity  model  (Figure  1);  b.  Red  line:  velocity 
estimate,  0%  data  noice  (Figure  7),  10  Gauss/Newton  steps;  c.  Grey  line:  velocity 
estimate  from  data  in  Fig.  10,  10  Gauss/Newton  steps. 


Figure  12:  Final  reflectivity  estimate  corresponding  to  Fig.  11. 
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Figure  13:  Gulf  of  Mexico  data:  CDP  stack. 


data  were  sorted  to  midpoint  order,  merged  into  bins  containing  up  to  301  offsets, 
interpolated  to  a  uniform  spacing  of  7.5  meters  and  radially  slant  stacked.  The 
antenna  pattern  of  the  receiver  array  was  computed  and  compensated  to  produce  a 
true  plane-wave  decomposition.  The  250  computed  plane  wave  slownesses  correspond 
to  an  angular  range  of  11  to  70  degrees  at  the  streamer. 

The  result  is  shown  in  Fig.  14.  Inspection  of  Fig.  14  shows  a  rich  sequence  of 
apparently  primary  reflections,  culminating  at  approximately  2.3  sec.  in  the  DHI 
event  mentioned  above.  Below  the  DHI,  the  character  of  the  data  changes  markedly, 
to  a  regime  of  slower  events,  presumably  multiple  reflections. 

An  additional  source  estimation  process  was  carried  out  after  the  plane-wave  de- 
compostion.  Both  wavelet  shape  and  directivity  were  re-estimated.  A  line  average 
of  the  water-bottom  reflection  event  was  computed  and  the  result  fit  to  the  wave 
equation  prediction  for  P-wave  reflection  at  a  fluid-solid  interface.  The  resulting 
p-dependent  estimated  source  wavelet  was  used  in  the  inversion. 

Finally  the  tau-p  gather  of  Fig.  14  was  subjected  to  a  deghosting  filter  and  a 
further  tapered  50  Hz.  high-cut  filter,  resampled  to  .004  sec.,  and  windowed  to  100 
traces  of  duration  2.5  sec.  The  resulting  gather,  shown  in  Fig.  15,  was  the  input  data 
set  for  the  inversion. 

Application  of  the  damped  Gauss-Newton  algorithm  produced  the  velocity  esti¬ 
mates  depicted  in  Fig.  16.  These  are  displayed  alongside  a  velocity  profile  extracted 
from  a  sonic  log  taken  from  a  well  near  the  midpoint  in  question. 

The  first  of  several  difficulties  in  assessing  these  results  is  evident  in  Fig.  16.  The 
estimated  velocity  compares  quite  well  with  the  log  velocity  trend,  where  both  are 
available.  Since  the  log  only  covers  a  segment  of  the  well  (in  this  case  approximately 
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Figure  14:  Gulf  of  Mexico  data:  Plane  wave  decomposition  of  a  midpoint  gather. 


TRACE  NUMBER 


Figure  15:  Gulf  of  Mexico  data:  windowed,  deghosted,  filtered,  and  resampled  version 
of  Figure  14.  Input  data  for  inversion. 
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Figure  16:  Results  of  inversion:  top  curves  are  velocity  estimate  (smooth  curve) 
and  interpreted  sonic  log  (ragged  curve).  Bottom  curve  is  stacked  reflection  series 
estimate.  Products  of  10  Gauss-Newton  steps. 

1400  -  2500  meters  depth),  however,  the  full  depth  range  of  the  velocity  estimate 
cannot  be  compared. 

Another  index  of  successful  velocity  estimation  is  the  extent  to  which  NMO  correc¬ 
tion  flattens  the  reflectivity  gather.  The  NMO-corrected  estimated  reflectivity  gather 
is  displayed  in  Fig.  17.  Several  features  of  this  gather  deserve  comment.  First,  the 
lateral  continuity  of  events  is  considerably  enhanced  relative  to  that  displayed  in  the 
data  gather  (Fig.  15).  This  moveout-  adapted  dip  filter  effect  was  noted  in  the  previ¬ 
ous  section.  Most  of  the  moveout  has  indeed  been  removed.  Some  residual  moveout 
remains,  however;  apparently  the  reflectivity  is  slightly  overcorrected  at  the  shallower 
depths  and  undercorrected  at  the  DHI.  We  speculate  that  this  effect  is  partly  due 
to  the  peculiar  structure  of  the  DHI:  careful  inspection  of  the  data  shows  that  this 
event  undergoes  a  phase  shift  at  about  trace  60,  and  another  at  about  trace  85.  We 
are  unsure  as  to  the  cause  of  these  phase  anomalies,  which  are  present  in  neighbor¬ 
ing  midpoints  as  well.  Obvious  possibilities  are  elastic  reflection  effects  and  thin-bed 
tuning.  In  any  case,  the  dip  filter  effect  of  our  optimization  turned  the  phase  shift 
into  moveout,  which  forced  overcorrection  of  overlying  events. 

Since  we  used  a  globally  precritical  gather  with  a  fixed  range  of  slowness  for  all 
times  (i.  e.  a  rectangular  mute  pattern)  the  range  of  angles  available  to  constrain 
the  shallow  paxt  of  the  velocity  is  quite  small.  We  expect  that  a  mute  pattern  better 
adapted  to  the  velocity  structure,  employing  different  slowness  ranges  at  different 
depths  to  keep  the  angular  aperture  as  constant  as  possible,  would  better  constrain 
the  shallower  velocities  and  prevent  this  overcompensation  for  deeper  non-primary- 
acoustic  artifacts,  to  some  extent.  Such  a  mute  pattern  could  be  estimated  adaptively; 
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Figure  17:  NMO  corrected  reflectivity  corresponding  to  results  of  Fig.  16. 

we  hope  to  implement  adaptive  mute  estimation  in  a  future  version  of  the  code. 

It  is  also  known  from  logging  that  the  density  fluctuates  considerably  throughout 
this  region.  The  general  acoustic  reflectivity  varies  systematically  with  slowness  in  a 
different  way  than  the  constant-density  reflectivity.  We  have  verified  in  subsequent 
numerical  experiments  that  neglect  of  this  difference  can  cause  systematic  misesti- 
mation  of  velocity,  through  the  amplitude  dependence  of  the  incoherency  operator 
W .  Since  our  current  code  is  based  on  the  constant-density  assumption,  some  of  the 
inaccuracy  of  the  velocity  estimate  might  be  due  to  actual  density  fluctuations  in  the 
real  earth.  This  is  an  obvious  matter  for  further  study. 

The  identities  and  depths  of  events  also  diagnose  the  success  of  an  inversion.  In 
the  present  approach  these  are  estimated  by  the  stacked  reflection  series  estimate, 
which  is  also  displayed  alongside  the  log  in  Fig.  16  (bottom  curve).  The  DHI  is  easily 
identifiable,  as  is  the  event  near  1520  meters  Both  are  placed  about  80  meters  too 
deep,  consistent  with  the  slightly  high  velocity  estimate  in  the  logged  zone  and  the 
overcorrection  of  shallow  events  noted  above. 

In  summary,  this  preliminary  test  achieved  a  modest  degree  of  success:  a  reason¬ 
able  velocity  model,  accounting  for  the  bulk  of  the  moveout  in  the  data,  was  estimated 
entirely  automatically.  Several  causes  for  the  residual  inaccuracies,  consistent  with 
the  layered  medium  hypothesis,  suggested  themselves,  and  remedies  will  be  pursued 
in  subsequent  work. 
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6  Coherency  Optimization  for  Laterally 
Heterogeneous  Models 

The  essential  ingredients  of  the  approach  to  velocity  inversion  sketched  in  the  previous 
section  were: 

(i)  parameterization  of  the  reflectivities  as  time  “sections”  (i.e.  traces),  so  that 
the  seismogram  is  a  regular  function  of  the  reflectivities,  one  reflectivity  (trace) 
per  plane- wave  “shot”  (synthetic  source); 

(ii)  referral  of  each  time-section  reflectivity  to  depth,  and  assessment  of  the  de¬ 
pendence  of  the  resulting  suite  of  depth  sections  on  the  “shot”  parameter  (i.e. 
slowness). 

In  this  section,  we  maintain  the  fiction  that  the  seismogram  is  adequately  ap¬ 
proximated  by  the  multi-dimensional  version  of  the  convolutional  model,  i.e.,  the  lin¬ 
earization  about  a  smooth  background  velocity,  and  moreover  apply  high-frequency 
asymptotics  freely.  A  natural  interpretation  of  (i)  and  (ii)  emerges  in  this  context. 

First  recall  why  the  parameterization  by  two-way  time  resulted  in  regularity  of 
the  seismogram  as  a  function  of  both  the  background  velocity  and  the  reflectivity. 
In  effect,  each  reflector  was  associated  with  the  time  of  arrival,  at  the  surface,  of  its 
reflection.  If  the  background  velocity  is  changed,  then  depth-parameterized  reflectors 
remain  fixed  while  their  reflection  times  change,  whereas  time-parameterized  reflectors 
remain  fixed  while  their  depths  change.  In  the  former  case,  high-frequency  arrivals  are 
time-shifted,  while  in  the  latter  case  they  are  not.  The  time-shift  is  a  time  derivative 
in  the  infinitesimal  limit,  and  its  appearance  marks  the  loss  of  regularity  of  the  depth- 
parameterized  reflectivity-to-travel-time  map.  For  time-parameterized  reflectivities, 
no  such  time  shift  occurs  as  a  result  of  background  velocity  change,  and  regularity  is 
maintained. 

Reflectors  and  reflection  arrivals  are  both  (near)-singularities,  i.e.  locii  of  high- 
frequency  energy.  In  the  high-frequency  limit,  therefore,  we  must  ask:  how  can 
we  parameterize  reflectivity  so  that  a  singularity  in  the  reparameterized  reflectivity 
corresponds,  under  conversion  to  ordinary  spatial  coordinates  and  mapping  to  the 
seismogram,  to  a  singularity  in  the  same  location?  In  several-dimensional  problems, 
singularities  may  have  orientations,  and  these  must  be  preserved  as  well. 

This  question  is  answered  by  a  theorem  of  Rakesh  (Symes  (1985),  Rakesh  (1988)), 
together  with  a  construction  presented  for  this  problem  by  Beylkin  (1985).  To  fix 
ideas,  suppose  we  consider  the  linearized  acoustic  problem  for  a  point  source,  which 
models  a  shot-gather: 

1  d2u  2  _  2 rz  d2uo 

U~~c[  ~dF 

where  rz  =  Cr/c,  is  the  “reflectivity  depth  section”,  cs  is  the  background  velocity,  u 
is  the  scattered  field,  and  u0  is  the  reference  field  satisfying 

^  “  V2 uo  =  /(05U  -  £.)> 
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2^  being  the  source  location.  Rakesh  showed  that  for  the  impulsive  case  ( f{t )  =  ${t)) 
a  singularity  in  rz  at  the  subsurface  location  y,  across  a  surface  element  with  normal 
2,  corresponds  to  a  singularity  in  the  reflected  field  at  receiver  location  xr,  time  tr, 
only  if  there  exist 

—  an  incident  ray  7 ,•  associated  with  the  reference  field,  emanating  from  the 
source- point  Xj  at  t  =  0; 

-  a  reflected  ray  7r,  passing  over  the  receiver  point  x^.  at  time  tT 

so  that  7,  and  7r  meet  at  the  reflector  point  y  at  some  intermediate  time,  making 
equal-angles  with  the  reflector  normal  rj_.  For  non-impulsive  but  broad-band  time 
signatures  f(t),  this  rule  governs  the  arrival  of  (primary-reflected)  high-frequency 
energy. 

This  is  exactly  the  usual  picture  of  the  reflection  process,  of  course.  Rakesh’s 
theorem,  which  justifies  this  picture  in  terms  of  the  wave  equation,  holds  in  complete 
generality,  so  long  as  the  background  velocity  c  is  smooth.  It  tells  us  exactly  which 
singularity  in  the  time  section  must  be  mapped  to  a  given  singularity  in  the  depth 
section  if  the  seismogram  is  to  return  such  singularities  to  their  original  position  and 
orientation.  Beylkin’s  construction,  on  the  other  hand,  applies  only  when  no  caustics 
exist  in  the  incident  field,  i.e.  each  depth  point  £  is  joined  to  the  source  point  x,  by  a 
unique  incident  ray.  Then  any  mapping  having  the  singularity-moving  properties  just 
described  must  differ  only  by  a  (possibly  frequency-dependent)  amplitude  modulation 
(technically,  a  pseudodifferential  operator)  from  the  Kirchhoff  migration  formula 

r2(y,xj  =  /<>(£, (g) 
=  J  dxr  tc(xJ,y,xl.)r(xr,r(x,,£,xr),xi). 

Here  t(xs,  j/,xr)  is  the  two-way  reflection  phase,  i.e.  the  time  from  x,  to  £  to  xr,  and 
w(xr,y,x,)  is  a  slowly- varying  amplitude  modulation. 

Thus:  the  reflectivity  time-sections  r(xr,<,x,)  will  be  converted  to  reflectivity 
depth  sections  r*(y,Xj)  via  a  Kirchhoff  migration  formula  like  (8).  This  relation 
implicitly  defines  the  seismogram  as  a  function  of  r(xr,f,x»),  and  guarantees  that  it 
is  regular  as  a  function  of  c4,  r.  Note  the  apparent  similarity  of  (8)  to  the  travel-time 
change-of- variables,  appropriate  in  the  layered  case  (formula  (A.l),  for  instance). 

The  second  ingredient  (ii)  in  the  coherency  optimization  approach  is  evidently  the 
condition  that  the  depth-parameterized  reflectivities  r*(£,xj  are  actually  indepen¬ 
dent  of  Xj  (“Every  shot  sees  the  same  earth”).  Regarding  the  source  locations  as  filling 
up  a  continuum,  this  amounts  to  the  condition  V£<r2(y,x,)  =  0.  The  composition 
rules  for  oscillatory  integrals  (the  local  calculus  of  Fourier  Integral  operators — e.g. 
Duistermaat  (1973),  Ch.  2)  give  the  result: 

V£jr,  =  V£<A>  =  I<(Vr,  +  P)r 

up  to  an  error  rapidly  decaying  in  frequency  content,  where  P  is  a  pseudodifferential 
operator  in  xr  and  t ,  i.e.  an  oscillatory  integral  of  the  form 

Pr(xT,t)  =  J  J  dxdsdkduj  eiM-r~-)+w{t~s)]p{xT,t,xs,k,uj)r(x,s). 
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The  amplitude  p  depends  on  the  ray  geometry,  i.e.  on  the  background  velocity  c4, 
but  the  phase  k(xT  -  x)  +  u(t  -  3)  does  not.  Thus  ( +  P)r  is  regular  as  a  function 
of  Cj,  i.e.  perturbing  cs  does  not  result  in  the  appearance  of  higher  derivatives  of  r, 
just  as  was  the  case  for  the  formula  (A.l)  for  W[c,,r]. 

Now  let  L  be  any  operator  from  depth-parameterized  reflectivities  to  sections 
which  has  the  reverse  effect  to  K  on  singularities:  for  instance  L  might  be  taken 
as  the  linearized  seismogram  operator  itself.  Then  another  application  of  the  same 
reasoning  shows  that 

W[ca,r]:=LV,Kr 

is  regular  as  a  function  of  c4,  r  for  the  same  reason.  We  take  this  formula  as  our  defini¬ 
tion  of  the  incoherency  for  the  laterally  heterogeneous  acoustic  problem.  Note  that  L 
must  be  computed  to  form  the  seismogram,  and  K  is  the  Kirchhoff  migration  operator 
(or  an  equivalent).  Thus  W[c4,r]  involves  only  well-understood  computations. 

The  attentive  reader  will  note  the  immediate  resemblance  to  the  definition  of 
W[c4,r]  given  in  Section  2. 

Thus  the  transformation  f  — ►  Kr  serves  the  role  of  the  travel-time  transforma¬ 
tion.  The  special  role  of  travel-time  in  regularizing  1-d  problems  was  noted  in  the 
work  of  Gray  (1980).  Some  time  later,  Gray  and  Hagin  (1984)  attempted  general¬ 
ized  travel-time  coordinates  for  laterally  heterogeneous  point  source  problems,  with 
limited  success;  travel-time  (ray-straightening)  coordinates  per  se  do  not  exist  in  gen¬ 
eral  for  several-dimensional  problems.  Nonetheless,  the  operator  K  accomplishes  the 
principal  effect  of  the  1-d  travel-time  coordinate,  i.e.  to  make  reflector  location  in¬ 
dependent  of  background  velocity  in  both  the  (reparameterized)  reflectivity  and  in 
the  seismogram  simultaneously.  Of  course,  K  is  a  more  complicated  operator  than  a 
change  of  coordinates,  except  for  1-d  (plane-wave,  layered  model)  problems. 

The  very  close  relation  of  the  method  described  here  to  “migration  velocity  anal¬ 
ysis”  (e.g.  Al-Yahya  (1989))  is  also  clear.  In  fact  W[c4,r]  vanishes  exactly  when 
the  migrated  gathers  Kr  are  flat  —  in  fact  all  traces  are  equal  —  with  respect  to 
shot  parameter,  i.e.  “when  sorted  to  receiver  gathers.”  The  differences  between  the 
two  approaches  are  exactly  parallel  to  the  differences  between  the  layered  coherency 
method  and  standard  NMO-velocity  analysis,  as  described  at  the  end  of  Section  3. 
Specifically,  the  method  described  here  is  based  on  a  differential  measure  of  coher¬ 
ence,  applied  to  the  model ,  whereas  migration  velocity  analysis  is  based  on  an  integral 
measure  (analogous  to  semblance)  applied  to  (migrated)  data. 

We  conclude  that  both  main  ingredients  of  coherency  optimization,  as  explained 
in  Section  2  for  layered  acoustics,  generalize  in  an  acceptable  way  to  a  simple  laterally 
heterogeneous  model.  Many  details  remain  to  be  settled,  some  of  which  will  doubtless 
be  crucial  to  computational  efficiency.  Also,  the  analysis  cited  in  Section  3  remains  to 
be  generalized.  A  technical  complication  is  that,  while  the  operators  L  and  K  exist 
in  general,  the  composition  rules  leading  to  the  regularity  of  the  incoherency  W[c„r ] 
must  be  modified  when  caustics  are  present.  Nonetheless,  we  have  established  that 
the  coherency  approach  is  not  restricted  conceptually  to  layered  problems. 
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T  Discussion  and  Conclusion 

7.1  The  scope  of  the  coherency  approach 

From  a  theoretical  point  of  view,  and  perhaps  from  a  practical  one  as  well,  the  chief 
defect  in  the  results  of  the  preceding  sections  on  the  layered  acoustic  problem  is  the 
absolute  lack  of  any  provision  for  multiple  reflections.  This  defect  is  cured  in  the  com¬ 
panion  paper  Symes  [1988c],  in  which  the  fully  nonlinear  bandlimited  layered  velocity 
inversion  problem  is  treated  with  full  mathematical  rigor.  We  reach  a  qualitatively 
identical  conclusion  about  the  appropriate  version  of  coherency  optimization:  that 
is,  it  gives  a  (regularized)  solution  of  the  inverse  problem  stably  dependent  on  near- 
consistent  data,  provided  that  sufficiently  many  reflectors  are  present,  i.e.  that  the 
target  profile  is  sufficiently  rough.  We  give  a  precise  sense  for  “rough”,  and  a  relation 
emerges  between  stability,  aperture,  bandlimits,  and  roughness  (reflector  density) 
very  similar  to  that  explained  in  Section  3. 

Both  convolutional  model  and  fully  nonlinear  versions  of  other  layered  inverse 
problems  should  succumb  to  the  same  approach.  We  mention  specifically  nonconstant- 
density  acoustics,  the  “marine”  elastic  model  (with  sources  and  receivers  in  an  over- 
lying  fluid  layer),  and  either  of  these  with  the  source  time-dependence  and  directivity 
also  regarded  as  unknowns.  Some  idea  of  the  novel  features  of  these  problems  in  con¬ 
volutional  approximation  may  be  gleaned  from  Sacks  and  Symes  (1987)  and  Bube, 
Lailly,  Sacks,  Santosa,  and  Symes  (1989).  No  technical  obstacles  appear  to  lie  in  the 
way  of  an  analogous  treatment  of  these  problems. 

Note  that  the  density,  regarded  as  independent  of  velocities,  will  not  be  recovered 
with  trend,  in  any  of  these  problems,  as  density  trend  perturbations  do  not  affect 
ray  geometry.  This  gross  density  ambiguity  is  well-known  (Tarantola  (1986))  and  has 
been  observed  in  output-least-squares  results  (Canadas  and  Kolb,  (1986)). 

The  program  outlined  in  Section  6  appears  feasible.  On  the  other  hand,  while 
an  analogous  coherency  optimization  principle  can  be  formulated  for  fully  nonlinear 
laterally  heterogeneous  models,  its  analysis  will  require  fundamental  advances  in  the 
understanding  of  wave  propagation  in  rough  media. 


7.2  Conclusion 

We  have  presented  a  novel  approach  to  the  reflection  seismic  inverse  problem,  which 
has  its  roots  in  utterly  commonplace  concepts  in  seismic  data  processing.  We  have 
formulated  this  coherency  optimization  principle  precisely  for  the  convolutional  ap¬ 
proximation  to  the  layered  constant-density  acoustic  model.  We  displayed  the  results 
of  numerical  experiments  with  both  simulated  and  field  plane-wave  data  sets.  The 
results  of  these  experiments  suggest  that  the  method  is  “practical,”  and  that  it  pro¬ 
duces  reasonably  accurate  and  stable  estimates  of  both  velocity  trends  and  reflectiv¬ 
ities,  to  the  extent  that  these  are  determined  by  the  precritical  plane-wave  data  set 
used.  In  particular,  the  coherency  method  extracts  reasonable  velocity  estimates  in 
cases  in  which  output  least  squares  inversion  fails  completely,  viz.  when  the  trend 
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of  the  initial  velocity  estimate  is  grossly  incorrect.  Finally,  we  formulated  an  anal¬ 
ogous  principle  for  laterally  heterogeneous  velocity  inversion,  a  problem  for  which 
output-least-squares  inversion  is  so  inefficient  as  to  be  infeasible.  It  remains  to  be 
seen  whether  coherency  optimization  yields  a  computationally  tractable  approach  to 
such  several-dimensional  problems. 
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